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Abstract 

This paper introduces the improved edge-based smoothed finite element method 
(bES-FEM) on simplices to solve the nearly incompressible elasticity problem. The 
improved method uses the piecewise linear polynomial space enriched with a bub- 
ble function on each element for the displacement, and the discontinuous pressure 
which is defined on a third mesh satisfies the uniform inf-suf condition. Moreover, 
the bES-FEM method induces a further softening to the bilinear form allowing 
the weakened weak (W2) procedure to get a quality solution. Several numerical 
examples are provided to show the effectiveness and reliability of the bES-FEM 
method. 

1 Introduction 

Elasticity materials such as rubber and rubber-like materials are used so much in 
the industry, we then have to research their structure and resort to approximate 
structural-analysis techniques. One of important features in the stress analysis of 
these materials is the nearly-incompressible case, i.e their bulk moduli are much 



larger than their shear moduli, or the Poisson ratio is near to one half. Computa- 
tionally, if we apply low-order finite elements based on quadrilaterals, hexahedra 
or simplices to solve a nearly incompressible elasticity problems, these methods 
are well-known to be the locking-effect. Until now, there are several numerical 
methods suggested to overcome this poor phenomena, for examples: the h-version 
finite elements [3111] result in poor observed convergence rates in the displacement, 
the B-bar method [15] for linear and nonlinear elasticity is leaded from separation 
of volumetric and deviatoric, the mixed formulations [H [8] having the inde- 
pendent displacement and pressure fields, enhanced assumed strain (EAS) modes 
[3U1 113 HI]) reduced integration stabilization [32] -[IS], two-field mixed stress el- 
ements [17] . Besides, we have several publications focusing on the analysis of 
average nodal pressure formulation based on enforcing a constant pressure field 
over a patch of triangles or tetrahedra [1H1 H9J EH EH ES]- It is easy to implement 
such scheme in the existing codes and has many nice properties. Although there 
exist many approaches to solve a nearly incompressible elasticity problem on the 
triangulation, a few methods are based on the rigorous mathematical analysis. 
For example [26] . the author introduced the method whose the pressure is discon- 
tinuous, and the displacement space belongs to the piecewise linear polynomial 
space of FEM enriched with a bubble function per element for the displacement. 
This method is based on the rigorous mathematical analysis. However, this dis- 
placement space certain has drawbacks of the FEM method: 1) overestimation 
of stiffness matrix especially for nearly incompressible and bending dominated 
problems, 2) poor performances when meshes are distorted, 3) poor accuracy for 
stresses. 

In this paper, we propose a new finite element method called by the edge-based 
smooth finite element method enriched by bubble functions (bES-FEM) for nearly 
incompressible elasticity problems, where the usual two types consisting of the 
£th-power bubble functions with £ = 3(2D), 4(3D) and the hat bubble functions 
are used. The bES-FEM method is designed in order to inherit advantages of 
some classical methods. Firstly, we use mixed methods [U [5] to reformulate from 
the linear elasticity problem ([!]) into the mixed displacement-pressure problem 
([5]), which leads to yield good approximations to the pressure as well [3J. The 
pressure is piecewise constants on the third mesh, the degree of freedom associat- 
ing with the pressure variable can be computed by linear combinations depending 
on the displacement variable. Secondly, the approximation displacement is com- 
bined by the edge-based smooth finite element method (ES-FEM) [23J and the 
bubble functions [28], [25]. The ES-FEM method is a member of the smooth finite 
element S-FEM family including the cell-based smoothed finite element method 
(CS-FEM) [21], the node-based smooth finite element method (NS-FEM) [22], the 
edge-based smooth finite element method (ES-FEM) and the face-based smooth 
finite element (FS-FEM) [21]. Moreover, in [9], the authors investigated the ex- 
tension of the ES-FEM method to 3D. The essential idea in the S-FEM family is 
to reconstruct the compatible strain field in finite element settings using the strain 
smoothing technique [12]. The reconstructed strain field is obtained over various 
sub-domains called smoothing domains created on cells (CS-FEM), nodes (NS- 
FEM), edges (ES-FEM) or faces (FS-FEM) of the mesh background, and the art 
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of the S-FEM family is the innovative design of the smoothing domains for desired 
amount of softening effects. The numerical operations used in the S-FEM family 
bring in information from the neighboring elements in desired ways. Depends on 
the requirements of the analyst, the property of an S-FEM model can be distinct 
in various ways from that of the standard FEM model. The S-FEM family can 
also be viewed as a combination of the numerical treatments of both used in the 
FEM and meshfree methods |19] . Since smoothing domains in the S-FEM family 
often (not always) cover part of the adjacent elements, the number of supporting 
nodes associated with the smoothing domain is larger than the number of nodes 
of an element. As a result, the bandwidth of the stiffness matrix (if formed) in 
the S-FEM family is increased, and the computational cost is hence higher than 
the standard FEM with the same sets of nodes. On the other hand, thanks to 
the propagation of non-local information brought by the adjacent elements, the 
S-FEM family often produces much more accurate solutions than the standard 
FEM method. For a given computational cost, all S-FEM models have better 
accuracy than that of the displacement-based FEM |15l I29j . The S-FEM models 
were applied to a wide range of practical mechanics problems [30 -[37J, and the 
S-FEM family has become a simple and effective tool for the analysis of a variety 
of practical problems. Among the existing S-FEM models, the ES-FEM was found 
so far the most computationally efficient [20J. Hence, the ES-FEM exhibits some 
interesting properties for solid mechanics problems such as: 1) it produces much 
more accurate solutions than the linear triangular elements (FEM-T3) and often 
found even more accurate than the FEM using quadrilateral elements (FEM-Q4) 
using the same sets of nodes; 2) the ES-FEM performs well with distorted meshes; 
3) it is super-accurate and has super-convergent properties in stress solutions; 4) 
the ES-FEM is stable even for dynamic analysis and 5) it is simple to implement 
into the existing FEM packages without any additional degrees of freedom. How- 
ever, if the displacement is only approximated by the ES-FEM method without 
enriching the bubbles functions, this method also is violated the in-suf condition, 
which is proved in this paper. Thirdly, using the bES-FEM method, we show that 
the inf-suf condition satisfies uniformly for the discontinuous pressure defined on 
a third mesh. The degree of freedom associating with the pressure variable can be 
computed by linear combinations depending on the displacement variable, which 
is a important feature be different from the classical mini element [5]. 

The organization of the rest of this paper is outlined as follows. In the next sec- 
tion, we briefly recall the boundary-value of linear elasticity problems, the mixed 
displacement-pressure form and weak form. Section 3 represents the description 
of ES-FEM method with bubble enrichment. Section 4 shows the mathemati- 
cal properties of the bES-FEM method. Displacement, energy and pressure error 
norms are defined in Section 5 for precise qualitative examination of various mod- 
els. Several numerical examples are presented in section 6. Section 7 concludes 
with some main remarks and discussions on directions for future work. 
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2 The boundary- value of linear elasticity 

We consider a static linear elasticity problem governed by equilibrium equations 
in a bounded domain SI C , rf = {2, 3} with the Lipschitz boundary d^l. 

— div a = / in 0,. (1) 

Displacements u are prescribed on the boundary d£l 

u = on T. (2) 

In the elasticity problem, the relations between the displacement field u, the strain 
field e and the stress field a are: 

• The compatibility relation between the strains and displacements 

= 1, d : Eij(u) = - (djUi + diUj) in (3) 

where di = ^- with (xi, ...,x d ) G R d , e(u) = [f%-(X)] i j=T ^. 

• The constitutive relations for elastic solids can be given as 

(Jij(u) = X5ij£kk(u) + 2fJ-£ij(u) in Q or a(u) = De(u) in Q, (4) 

where /i = ' ^ = ( i+J)fi-2v) are * ne Lame constants and 5ij is the 

Kronecker delta index and D is the Hooke elastic matrix. 

Our attention is devoted to the nearly incompressible case which corresponds to v 
being very near to 0.5. This case leads to the poor performance that is well known 
to be the looking effect and instability. 

2.1 Mixed displacement-pressure form and weak form 

The elasticity problem ([!]) can be rewritten by the following mixed displacement- 
pressure form 

— div a = f in Q, (5) 

div u — — = in 17, (6) 
A 

where a pressure p is introduced as a extra variable. The mixed form is the same 
as penalized Stokes equations. 

The mixed approach finds a displacement field u G V = [M 1 and a pressure 
p £ Ll{Q) := G L 2 (Vt) : J q dx = oj that satisfy 

a(u,v) + b(v,p) = (f,v) VveVo, (7) 
b(q,u)-jc(p,q) = VgGL^O), (8) 



4 



where p belongs to Lq(Q) because of the Dirichlet boundary condition, and the 
bilinear forms are defined, as follows: 

a(u,v) = 2fi J e T (u(x))De(v(x)) dx, 
h 

&(*,«) = j q{ x)V.u { x)dx, 
n 

c(q,p) = J q(x)p(x) dx, 
h 

(f,v) = / f(x)v(x) dx, 



with / G L 2 (J7), the matrix D of material constants is symmetric, positive definite 
and its eigenvalues are bounded in [a,b] C M + . 



3 Description of ES-FEM method with bub- 
ble enrichment 

3.1 The meshes and the finite spaces 

The polygon domain is discretized by the triangulation Th (the primal mesh), 
where Th consists of simplices, either triangles or tetrahedra. The set Th has N e 

_ N e _ 

elements, N n nodes (or vertices), N s edges and f2 = |J Tj. For each element 

i=l,TiET h 

T £ Th, the barycentric point ct is called a mesh point of T. Let Vh be the 
standard linear finite element space defined on the triangulation Th, 



V h = {ue [H\n)] d , u\ T e [P 1 (n)] d ] 



which has the standard nodal basic functions iVj (i = l,N n ) associated with the 
node i. 

Defining the space of bubble functions as 



Bh 



b T e H 1 ^), b T \ dT = and J b T (x)dx > 0,T G T h 



T 



where the basis bubble functions are considered in two types (see [25] and 
on each element T S Th, 

• The first one is the £th-power bubble function with £ = d + 1 

{d+i 
c b {d+l) 3 n A TW (x), if a; G T C T h , ^ 
elsewhere 

where each function X T (i) is a barycentric coordinate associated with a vertex 
x T (i) of the triangle T, and c& is computed in such a way that 6t(ct) = 1 
with the centroid ct of T. 
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• The second one is the hat function on T, where T is partitioned into sub- 
triangles (2D) or sub-tetrahedra (3D) {T^} i= l d+l by joining the centroid 
ct to the two vertices (2D) or the three vertices (3D) 



b T {x) 



c b (d+ 1)At w (x), if x € T (i) C T C 
elsewhere 



(10) 



Using the bubble space Bh, the finite element space for the displacement which is 
enriched with bubble functions is defined as Vjf = Vh [B h ] d C [H 1 (0)] <i . Each 
function 6 V® which is restricted on T 6 Th is written as 



u h (x) 



d+l 

£ 



diag N T (i)(x), N T (i)(x) u T{l) +diag 



\ 



N b c Jx),...,N b c Jx) 



d dimensions 



V 



II 



d dimensions 



/ 



e h (x)ev h 



b h (x)e[B h ] d 



(11) 

where the diagonal matrix is denoted by diag(.), for each i = l,d+ 1, iV^i) is 
the standard nodal basic function associated with the vertex x T (i) of the triangle 
T, N b T is the standard nodal basic cubic bubble function defined on T with the 
centroid ct- Besides, the values u T (0 and u CT G R d are the nodal values of at 
the vertex and the barycenter ct- 



Furthermore, in order to design a smoothed mesh for smoothing strain and 
divergence operators, this smoothed mesh that is called a dual mesh T h * is con- 
structed by connecting all vertices and mesh points of Th- The dual mesh T h * 



satisfies £1 



(J Of , and all elements of T h * are not overlapping together. 



fc=i,njerjj 



Each element Qf, G T h * is associated with an edge of the primal mesh Th ■ For 
examples, we have an illustration of an element corresponding to an inner edge 
efc 





Li: >--\7 


\ □ (••)•' / : ; ^ 

\ \ K" S 2/\ '. 


\ n / 










6] is a boundary edge 


on r 


e 2 'is an internal edge 





Qj,Q, are two elements of the dual mesh 

Figure 1 




A center of a face 

center of a tetrahedron 



and another illustration of an element Of located at the boundary 
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2 



Figure 2 



With the dual mesh 7^*, the space V® is equipped with the following the inner 
product, semi- norm and norm (see [T7]): 



(W, V) V B 



w\ 



\W\\ V B 



n Vl=1 7 



fe=i 



d d-l d 

E 0« (wi)9u ( v i) + E E (&j ( w < ) + (&i ^ ) + 5ji ( «j ) ) 

1=1 jr'=j+l 



fc=l 



1=1 

d 



d-l d 

^liiO^O+J^ Y i9ij( w i) +9ji(wj)Y 

i=l j=i+l 



i=l 



, / d \ Ar s 
/ \Y W i) dx + Y meas ( n k) 



d d-l d 

Y (9ij( W i) +9ji(Wj)Y 
i=l j=i+l 



i=l 



where vectors to = (toi, wa), v = (vx,...,Vd) belong to Vjjf, ft 
meaI(o s ) / ^ i ( x ) c ^ x on a smooth element Q|. £ 7^*, and the measure meas(0|) is 

equal to J ldx. 

Thank to the remarks 3.4 and 3.5 of [T7], we have relationships between |.|ys, 
||.|| v b and |.|i, ||.||i, as follows: 



\w\ v b < \w\\ and |M| v b < IMIi with w G C H 1 ^) 
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where which is a Sobolev space has the semi-norm |.|i and the norm ||.||i 

d 

defined from the inner product (u;, v)\ = ^ (u^,i>i)i (see chapter 3 in [27]). 

i=l 

Next, a third mesh 7^** is constructed by connecting all bary centric points 
{ct}t£Th an d midpoints of all edges of Th in 2D, plus barycentric points of all 

_ N n _ 

faces in 3D. The third mesh T^* satisfies O = |J Vi, and all elements of 



7 



7^* are not overlapping together. Each element Vk G 7^** is also associated with a 
vertex of the primal mesh. 






» k 












e i 






/ X ^ 1 









The polygon (x k ,X e ,X f ,X e ,X A ,X^,X f : ,C T ^ 
Thepolygon (c r , .v r , c ; . , ^ .c 3 . , jc^, , c 3 . , ) = ^ X , JC , X are midpoints of edges 

The polygon (c^.J^,^,^) = n2] X^,X^,X^ are barycentric points on faces 
Figure 3 

Basing on this third mesh, we define the following finite element space for the 
pressure 

V** = {pe L 2 (Q) such thatp \ V €V°(V) ,Ve T%* } , 

N n 

if ph G V", the function = ]C P«X«) where Xi are the characteristic functions of 

i=l 

Vi G 7^**, i = l,N n , and the norm ||.||o of V" is defined \\q\\o = (j q 2 dx^j~ with 

^ g y** ^ Now, we represent the bES-FEM method for the nearly incompressible 
elasticity problem: 

3.2 Smooth strain and smooth divergence 

On each smooth element Qf. G T£, the strain e(uh) is smoothed as 



1 



meas{Vt s k ) 



e(uh)(x)dx with G Vjf, 



(13) 



and for the smoothed divergence, one has 

(V.ith) |n? = / V.Uh{x)dx with u h eV®. 

' k meas(0, s k ) J 



(14) 



3.3 Weakened weak statement for the bES-FEM 

Here, we want to find the discrete solution (uh,Ph) G Vff x V£* such that 



a(u h ,v h ) + b(v h ,p h ) = (f,v h ) with Vv h G Vf , 
H v h,Ph) ~ {c(p h , q h ) = with \% G V ( 



(b) 



(a) 



(15) 



N s _ 

where a(u h ,v h ) = 2/i ^ meas(Q 8 k )e^ k ' (u h )D^ k \v h ), b(v h ,p h ) = J (V.u h )p h dx, 

k=i n 
c(ph,Qh) = fPhQhdx and (f,v h ) = f f(x)v h (x)dx. 
n n 

4 The mathematical properties 

Now, we show the following important properties of the bES-FEM method, when 
this method is applied into the linear elasticity problem: 

The first property: The bilinear form a(., .) is continuous, symmetric and semi- 
elliptic on Vf := {v G V® : b(v,q) =0,g6V*C i.e. there exists an 
«o > such that 

a(v,v) > a \v\ V B, veV® . 
Thanks to the theorem 3.1 in [18J, it help us to prove this property. 

The second property: The bilinear form &(., .) on V® x V^* is continuous and 
satisfies the uniform inf-suf condition, i.e there exists a positive constant /3q inde- 
pendent the mesh size such that 

Kuh,qh) . b{u h ,q h ) 

sup -7, — n > sup — | — r— > /Solkhllo, q h £V h . (16) 

Uh&V B \\U h \\ v B because of jl2} Ufi&V B \\U h \\l 

Proof 

Let (u h ,q h ) G x V£*, we have 

b(uh, qh) = j (V .u h )q h (x)dx = j V.(4 + b h ) q h (x)dx = J (V.4 + V.b h ) q h (x)dx, 
n n n 

(17) 

where there exist uniquely 4 G 14 and bh G such that u/j = 4 + bh- In 



(17), the smoothed divergences V.4 and which are restricted on Q s k G T h * 
are defined 

(VA) |ns = /VA(^ and (VA) |,, = /W)efc. 



To prove b(uh,qh) satisfying the uniform inf-suf condition, we need to look for a 

relationship between b(uh,qh) = J V .Uh(x)qh(x)dx and b(uh,qh) with (uh,qh) £ 

ft 

x V?*, because in [26] the uniform inf-suf condition holds for b(uh,qh)- 



In the first step, we compute the value of j (V .£h)qh(x)dx — j (V .£h)qh(x)dx with 

n n _ 

(4) QTi) S VfcX V£*. Using the fact that V.4 is constant on each T G 7/j, we obtain 

j Vl h q h (x)dx= (V-4) It / ^0*0^ (18) 
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At any element T G Th with its vertices {x T a)} i= 1 d+1 , we have 



/d+1 d+1 
q h {x)dx = ^meas (v Xt( . } D T) .q T{l) = (V-4) \ T ^ 
, i=l i=l 



meas(T).card(£ T (i) 



2(d+l) 

(19) 

where for each i = l,d + l, Kc (i) £ Tj** is respectively associated with a ver- 
tex x T (i) of T, and g T (0 is a nodal value of qn at a vertex x T (i). The value 



meas(K Cr( . ) (IT) is equal to 



meas(T).card(£ rr ^ ) 



2 (d + i) with i = 1 , d + 1 , because of the 

third mesh T£* constructed by barycentric points of all faces, midpoints of all 
edges and the centroid points ct for all T G Th- The set £ T ^ contains all edges of 
T such that these edges have a common vertex x T (i) . Additionally, for all i = 1, d, 
the notation card(£ T (i) ) which is the number of all elements of £ T (i> is equal to d, 
because the primal mesh Th is a triangulation. 



In (17), we consider 



(V.£ h )q h (x)dx = I (Vl h ).q h (x)dx 
T&T h J T 



(20) 



-q T (i), 



On the above element T E Th, the integral J (V .£h)-qh{x)dx is computed by 

T 



,. d+1 

/ (ylh).q h {x)dx = 

rr 1 = 1 



E 

e T(i) e£ T(l) 



meas [V x ,., nTnfZ; 



T (i) 



meas I 



Vl h dx 



r(i) 



■ q T (i). 

'21) 



where the domain VL s e £ T^ corresponds to the edge e T (i) . 



In the first case of T (a triangle or tetrahedron), we assume that all edges of 
T are inner edges, i.e the edges are not on the boundary d£l. For each i = 1, d + 1 
and j = l,d, the integral J V.lhdx is computed by 

\(.) 

/" V.l h dx = meas n r) . (V-4) |t+ meas (fi^ n . (V-4) |k , 



Or 



Ke% T(i) \{T} 



(22) 



where the notation T e is a subset of 7h such that its elements have a common 
edge e r( i) and T G X ... . 



From (21 ) and (22 ), the integral J (V.£h)-Qh{x)dx has the coefficient of (V.4) It 

T 



E 

e T(i) ££ T(l) 



meas[V x ...nrnflf , J .meos fif , nT 



meas ( f2f , , 

e T( l ) 



(23) 
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Together J (V .£h)-Qh(x)dx , we only find the coefficient of (V.£h) \t QtW m / -^h)-Qh{x)dx 



for all K G T e (i) \{T} and e T (*) £ ^tWj as follows: 



meas ( V^ t(i) nATlOf ,., ) .meas ( ftf ,., fl T 



n(i) 



meas ( 



"(0 



(24) 

From (23 ) and (24), in the integral J (V .£h)qf l {x)dx, the coefficient of (V.£h) It Q'tW 

is equaTto 



E 



e T(i) e£ T(i) 



meas K ... nTnO| , % .meas i Qi , nT 



meas I f2f 



"(0 



+ 



meas Vi ,,n.ft"nrH ,. s .meas fi= , nT 



^^UT} 



meas Q?, 



(25) 



By using the centroids ct for all T E 7^, the midpoints of all edges, plus barycentric 
points of all faces (3D) to construct the dual mesh 7^* and the third mesh 7^**, we 
have 



= , , / , meas ( V x ... nKntt s p 
2(d+ 1) V x tM 



E 



meas(L) 



, „„ _\ meas(T) / , 

meas O nT = — , — — , raeos 11 , , 

6 tM / d+1 V e T«y ^ — < d+i 



Therefore, the coefficient of (V.%) |t 9tM ^ s computed by 



meas(T) 
d + 1 



E 



s T(i) e£ rW 



meas(T) 
2(d+l) 



\—\ meas(L) 

.E d+i 

LeT: 



E 

K€T« \{T) 



meas(K) 
2(d+l) 



E 



meas(L) 
d+1 



e t(0 



meas(K) 
2(d+l) ' 



(26) 



meas(T) .card(£ T (i 
2(d+l) 



(27) 



From (19) and (27), the two coefficients of (V.^) |t <2tW m the ^ wo integrals 
J (V .£h)Qh(%)dx and J" (V .£h)Qh{x)dx are together equal, 
n n 

In the other cases of T, it has least one edge belonging to the boundary dQ, we also 

obtain the same results. Hence, the value of J (V .£h)lh(x)dx — J (V .£h)Qh(x)dx 

n n 

is equal to 0. According to the result of [26j, without enriching the displacement 
space, the uniform inf-suf condition is violated. 
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In the second step, we find the relationship between J (V .bh)Qh(x)dx and J (V ' .bh)qh{x)dx. 

n n 
Using the definitions of the spaces Bh and V", with (bh, qh) £ Bh x V", we get 



I V .b h (x)q h (x)dx = ^ ^ 



X7.b h (x)qidx 



(28) 



and 



/ V .b h (x)q h (x)dx = ^ ^ 
j! TeT h i=i,v<eT£ 



Tnv^0 



V.b h (x)qidx 



WiDT 



(29) 



Considering T whose all edges stay in the internal domain f2, for each i = 1, d + 1, 
we have 



V ' .bh{x)q T (i)dx = q T( i)U CT . / VN^ T (x)dx, 
■/, is rcwnt icii as &/, = U Cj ,l. Cl 



(30) 



Vi m nT 



where 6/,, is rewritten as 6/,, = u r ^N„„ with a vector u CT , and 



V .bh{x)q T (i)dx 



E 

e w e£ w 



E 



meas K n nTnfi 
V >(«) 



s 



meas [0.1 , ^ 



E 



V.bh(x)dx 



02 ,.,r\K- 



meas [ O 



E 



/ 



VN*(x)dx 



<ZtW 

\ 
/ 



From (31), in the integral J (V .bh)q T (i)dx, the coefficient of q T (i)U CT with 



T G % T (ij is equal to 



E 



meas [ V x ... D T D 



r(i) 



meas ( 



fif flT 

TW 



(32) 



Furthermore, the other coefficients of qj>(i)U CT , which are also found in 



(y.b h )q T (i 



v x ..ok 
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are equal to 



meas [V x n K(l Q s p , , 



meas ( 



VNj!(x)dx 



(33) 



v - "' I jre{r e (1) \T}cr h 



Ve T(l) e£ T(i) 



Thank to ( 26 ) , ( 32 ) and ( 33 ) , the coefficient of q T (i) u CT are computed in the integral 



J (V .bf l )qf l (x)dx, as follows: 
n 



E 

c tW 6£ t(0 



meas f2S 



e T (i) 



+ 



E 

ft{ r v(.)\ T } c7; ' 



t( ! ) tW 



meas I Q 



e T (i) 



VN°(x)dx 



«S ,,nr 



E 

e T (i) e£ T(*) 



2(d+l) 



3+1 

TO) 



E 

«{ r v(.)\ T }^ 



/ 



meas(K) 
2(d+l) 



measfL) 



V ie7i rW 



ViV c 6 T (x)dx 



,,nT 



\ E 



VN b CT (x)dx 



\n% nT 



(34) 



In two dimensions, we compute the coefficient of q T {i)U CT (34) on the following 
triangle T having three vertices {xj, Xj, xj~} 





Figure 5 
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This coefficient is equal to 



J VN b CT (x)dx+ J VN b CT {x)dx =~ J 



\ I 
1 

2 



(if ,nT 



/ 



Nc T (x)n mdx+ I N b CT {x) njl) dx 



7?> 



(35) 




Figure 6 

In Figure 6, we introduce some extra notations including that the midpoints of 
edges [xj,Xj], [xk,Xi] and are denoted by x^, Xki and xjy. Besides, nota- 

tions j£\ 7^ 2) , jj 1 ^ and 7^ are edges [x k ,c T ], [xij,c T ], [xj,c T ] and [x k i,c T ]- Vec- 
tors n (i), n (2), n m and n r 2 ) are outward normal vectors of Of , DT, T4- flT, 

7fc 7fc 7} 7} Fi>%'J 1 

a,.] H T and n T, respectively. Length of each vector n (1), n (2), n (1) and 

(!) „,( 2 ) fl) _ , (2) 



n 



7 



(2) is equal length of each segment 7} , -f, , 7 ■ and 7 , so n (1) = 2n^ 



(2) 



and 



n (i) = 2ra (2), because length of each segment 7[ 1 ' ) , 7- ' is equal two times length 
of each segment 7^ , 7 • , respectively. 



7fc 



We directly compute the coefficient (35) in the following two types of bubble 
functions 

■ The ^th-power bubble functions M) with £ = 3 (the cubic bubble functions) 



Assume that T is the reference triangle, Mt is the Jacobian of transformation from 
the triangle T to f, J T = det(M T ), 9^ = f N b ,(x)dx, ef ] = f N b .(x)dx, 



6^ = J N b , (x) dx and 9^ = J N b . (x) dx, where notations jf* are defined 



(i) 



-(2) 
7 2 



4 1] 



-(2) 
73 



7i 



(i) 



r 1 -(2) _ r 1 -X 1 ) _ r 1 

[ X f (1) ' C f J' '1 — l X f (23) ' C fJ ' ^2 — l X f( 2 ) ' C TJ ' 



"(2) _ r I "(1) _ r I -(2) _ r i 

7 2 — [Xf( 13 ), Cf\, 7 3 — [Xf ( 3 ) , Cf\ , 7 3 — [Xf( 12 ),Cf\ 
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with points x^ (1) (0, 1), a^ {a) (0, 0) , x^ (3) (l,0), Xf (12) (0, |), Xf (23) (|, 0), Xy ( i 3 )(§, |) 

cUld C^p ( g ; 3 ) * 

^,(0,1) 




W0,0) 



^.3,(1,0) 



Figure 7 

Together this assumption, we use the lemma 3.2 of [22] to obtain 



N b CT (x) n 7 (D dx + / iV c b T (x) n 7 (D dx = J T (o^ h.p + 6^ ) M" 1 . (36) 



7« 



VN b c Jx)dx = J T [9 .(2)71 (2) + (2)77. 



72 7 2 



7 3 7 3 



(2) 



M 



(37) 



14, nT 



By directly computing the quantities on the reference element T, we have 

• The barycentric coordinates of a point P^ 1 ), x^ 2 )) in the reference triangle 
T are Ai(x) = A2(x) = 1 — x*- 1 ) — x^ 2 ) and \^{x) = xW with x = 
(x^ 1 -*, x^ 2 )). The basic cubic bubble function on the reference triangle T is 
iV c fe .(x) = 27.A 1 (x).A 2 (x).A 3 (x). 

• The segments Ji , 7i are on the line (di) x^ 2 ** = —2x^ + 1. 

• The segments 7 2 , 72 are 011 the line (g^) x^ = x^ 1 ). 

• The segments 73 are on the line (g?2) x^ = — 0.5x^ + 0.5. 

• The values of the following integrals: 

1/3 

\/5 



7i 



(l) 



.-,(2) 



7i 



7V c 6 T (x)d 7 (x) 



iV c fe T (x)d 7 (x) 



27 y s 2 (l - 2s)V5ds = 


1/2 

27 | s 2 (l- 2s h/5cfe = ^p 



9 - 16 f) 



-(2) 
7i 



1/3 



(38) 
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Mi) 

72 



7,(2) 
72 



1/3 

y JV*,(x)d7(a;) = 27 J s 2 (l-2s)V2ds 





1/2 



y JVj r (a;)d7(a?) = 27 J s 2 (l-2s)V2ds 



27y/2 
162 ' 

11.27.V2 
2592 ' 



9 - 16 a 



7 < a > 



1/3 



(39) 



73 



.,L, - / ^ T (x)d 7 (x) 



1 

27 y s(0.5-0.5s)Vl + 0.25ds 



16.VH25 

48 



73 



1/3 
1/3 



7s 



= / JV*,(x)d7(a;) 



27 y s(0.5-0.5s) 2 Vl + 0.25ds 



11.VT25 



-(2) 
7 3 



The relationships between the normal vectors ra„(i) and n„(a) with i = 1,3: 

7i 7j 

n ( i)=2n ( 2), n ( i)=2n ( 2) and n ( i)=2n ( 2). (41) 

7i 7j 72 7 2 73 73 



From (36)- (41), we point out that 



y VN b CT (x)dx+ J VN b CT (x)dx ' U 
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y VN b CT (x)dx. (42) 



fif ,nT 



Hence, we use the results of (30), (35) and (42) to imply that 



16 

f/ie coefficient of p T (i)U CT in b(uh,Ph) = —-the coefficient of p T a)U CT in b(uh,Ph)- 

(43) 

(44) 



With computations of ( 28 ) , ( 29 ) and ( 43 ) , we conclude that 

16 f 

(V .b h )q h (x)dx = — / (V .b h )q h {x)dx 



Defining = £h + j§bh, using (44) and the result of the first step, we get 



(V .u* h )q h (x)dx = / (V.u h )q h (x)dx. 



(45) 



Finally, thank to the result of Theorem 3.1 of and (45), the uniform inf-suf 
condition holds for the bilinear form b(., .) on x V^*. 



The hat bubble functions (10) 
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For each triangle T G Th, the divergence of hat bubble functions is equal to 
constant on each sub-triangle {Tu)}^ of T, so we get 

J VN b CT (x)dx+ J VN b CT (x)dx = ^ f VN b CT (x)dx. (46) 



By (30), (35) and (46), we obtain 



the coefficient of p T (i)U CT in b(uh,Ph) = the coefficient of p T ^)U CT in b(uh,Ph) 



which implies that 



(47) 



(V .b h )q h {x)dx = / (V .b h )q h (x)dx 



Therefore, 



(V .u h )q h (x)dx = I (V.u h )q h (x)dx. 



(48) 



(49) 



In three dimensions, we also compute the coefficient of q T a)U CT (34) on the follow 



ing tetrahedron T constructed from four vertices x T (2), x T (s), x T (4,)}, 

x 




Figure 8 

where ct is the centroid of T, x^P is the barycentric point of the triangular face 

(x T a) , x T (j) , x T (k) ) , Xe is the midpoint of the edge [x T (i),x T (j)], where i,j,k be- 
long to {1,2,3,4}. 
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In particular case, the detailed coefficient of q T (i)U CT (34) of b(u, q) is computed as 

/ \ 

J VN b CT {x)dx+ j VN b CT (x)dx+ J VN b CT {x)dx 



\ PtW^tW)] r T (i)' x T(fc)J Lvco'VwJ 



/ 



c T y f ',x tU) ) (c t ,z} > t(j -) 



(iji) \ V / tWJ ( (ikl) \ V / tWJ 



(50 



where j,k,l £ {1,2,3,4} are different from i. Vectors ni ( iik ) \i n ( WO 
ni Mia \, ni (iti) \ , n/ (y0 \, ni u kV) \ whose length is 

[CT,xy ,x T ( k )) { c Ty f ' x T (k)J [ c T,xy ' ,x t(1) ) [c T ,x K f ',x T (i)) 
equal to measure of triangular faces [c T ,xf k \x^)), (c T , xf \ x tU) ) , (c T ,xf k \x T{k) 

ct, Xj kl \ x T ( k )^j , \or,x^ l \x T (i)\ and (ct, Xj kl \ x T (i)^j respectively, are the out- 
ward normal vectors of T. 



Moreover, we also have the coefficient of q T (i)U CT of b(u,q), as follows: 



/ / NcT^nf^^^^^djix) + J N CT^ x ) n ( CT ^ l \ x (^)\ d ^^ x ) + \ 



(ijk) (i,jp 



c T ,xy ,x. 



CT,Xf ,X 



I N b CT {x)nf T ^ l]k) ^ k) ^d>y{x) + J N b T (x)n ( ^ ^ kl) ^ k) ^dj(x)+ 



c T ,xy ,x. 



(ikl) (i,k) 



ct,Xj : x 



Ct,X\ ,Xe 



.WO AW) 



c Ti x f i x 



(ikl) 



c Ti x f i x 



CT, x f i x e 



Furthermore, we have relationships between normal vectors in the two formulas 



(50) and (51) 



(51) 



[C T ,X) >,X t(J) 



2n (c T ,xf k \x^y n (c T> xf l \x T ^ 2 " 



CT^lx^' 



n ( (ijk) 

{c T y f ,x T ( k) 



In/ uj k ) u, k )\,n/ tiki) \ — In i u k \\ u, k )\, 
\c T ,xy ',xl' >\ \c T ,x^ f ',* T (fe) ) [c T ,x y f ',x),'' ] 



n ( WO \ 

[c T ,xy >,x T ( l} ) 



2n/ uji) (n)\,nf uki) \ —2n/ uki) un 
I c T ,xy ',xi ■ ') U T ,xy ',x T (i)\ ic T ,xy ' ,x). 



(52) 
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Obviously, from (50), (51) and (52), if the hat bubble function N^ T is considered, 



we also get 

the coefficient of ' p T (i)U CT in b(uh,ph) = the coefficient of p T (i)U OT in b(uh,Ph 
which leads to 



(V ' .b h )q h {x)dx = J (V.b h )q h (x)dx. 
n n 



(53) 
(54) 



Hence, 



(55) 



(V .u h )q h (x)dx = I (V.u h )q h (x)dx. 
a a 
With the 4th-bubble functions N^ T , there exists a positive constant a such that 



with Uh 



(V .u* h )q h {x)dx = J (V.u h )q h (x)dx. 
n a 

bh and u* h = + abh in $7. 



(56) 



The third property: The bilinear form c(.,.) is continuous, symmetric and 
positive semidefinite, i.e 

c(<Z,<z)>0, 9 G Lg(O). 

Remark |4j. 1: Prom (|15[)), we can compute the pressure field by the displacement 
based on the edge-based smoothing domains, as follows: 



J p h qhdx = X j (Vu h )q h dx 



For each i = i,N n , we choose q^ = %ii where Xi is a characteristic function of 



Vi 6 T£*. It follows that 



A 



Pi 



meas{Vi) 



meas(Vinn%)(y.u h ) |nj, 



(57) 



fe=i, fi£e7£ 



N„ 

E 

i=i 

transformed into 



where ph = YlPiXii Xi-Xj = with i ^ j. Then the bilinear b(vh,Ph) can be 



( 



b{v h ,p h ) 



" n i 
i — measiVi) 



^2 meas(Vi n£l s k ) (V.u)\ 



( 



Ns 



y~* j meas(Vi n Of 
=i,nfer* 

(58) 



V.v) 
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Therefore, we arrive at a problem of finding G such that 
a(u h ,v h ) + 



i 



E 



J meas(Vi) 
(/,0 V^eVf, 



V] meas(Vi n 0|) (V.it) 



fc=i,njeTjj 



\ 



^ meas(Vi n Of) (V.u)| n( 

,nfeT 



where the solution u/j of (59) is the same as the solution of the problem (15). 



Remark |4|2: For applying the bES-FEM scheme into the linear elasticity prob- 
lem, the algebraic formulation of the problem is written as 



A B T 
B -\C 



Uh 

Ph 



fh 




(60) 



where A, B, C are matrices associated with the bilinear form as a(., .), &(., .) 
and c(., .), respectively, and fh is associated with the linear operator (/, .). This 
framework of the bES-FEM scheme implementing in the linear elasticity problem 



is as in the case of the mini scheme. However, the matrix C of (60) is different from 



the associated matrix of the mini scheme, because the matrix C of (60) is diagonal 



and each degree of freedom corresponding to the pressure can be computed by 



(57). It follows that the matrix is positive definite. 



/ 



(59) 



5 Error norms 

Here, we introduce the notations of the well-known methods used to compare with 
the bES-FEM method— the edge based smooth finite element method enriched by 
bubble functions. The list of these notations includes 

• MINI - The mixed displacement /pressure finite element method with cubic 
bubble functions [2]. 

• FEM - The standard FEM using three node triangular element with shape 
linear function |29j . 

• NS-FEM - The node-based SFEM [22] using triangular elements. 

• ES-FEM - The edge-based SFEM |23j using triangular elements. 

• ASOI - Assumed strain stabilization using the 01 strain field |45j . 

• ASQBI - Assumed strain stabilization using the QBI strain filed [45] . 

• ASMD - Assumed strain stabilization using the strain field associated with 
the mean dilatation element |45j . 

• SRI - The selective-reduced integration element [TT] , 
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Furthermore, for studying the error and convergence rate of above numerical meth- 
ods, we offer three types of error norms consisting of the displacement error norm 
and the pressure norm. 

• Displacement error norm 

The displacement error norm is defined by 

1/2 



\ u - Uh\\L 2 (n) 



^ / {u - u h ) T (u - u h )dx 



(61) 



jeT h T 

where u is the analytical solution for the displacement and u/j is the approx- 
imate solution for the displacement of any numerical methods. 

• Pressure error norm 

We also have the following definition for the pressure error norm 

1/2 



\P~Ph\\L2(n) 



E /<p- 



Phfdx 



(62) 



where p, ph are the analytical pressure solution and the numerical pressure 
solution, respectively. 

• Energy error norm 

Recalling that we only use and approximate the displacement field in the NS- 
FEM method, but in the the MINI and bES-FEM methods, these method 
have not only the displacement variable but also the extra discrete pressure 
variable. Hence, we propose a different definition of the energy error norm 
for each method, as follows: 

For the NS-FEM, ES-FEM methods, the evaluation of the energy error norm 
is based on N s smoothing domain Of S T£ , so it is defined as in [18] 

■ 1/2 

\u- u/.Ul-^ rneas^l s k )[a-a^{u h )\ T D- l [a-a^\u h )] 

k=l 

J 

(63) 

where a is the analytical solution for the stresses and a^ k \uh) = De^ k \uh) 
is derived from smoothed strain solutions on smoothing domains fit. 



For the bES-FEM method, it is also based on N s smoothing domain G 7^*, 
but the definition of its energy error norm depends on the displacement and 
pressure 

Ns "v V2 

2/i £ meas^l)[e(u)-^ k \u h )} T D^[e{u)-e^{u h )] 



\U-U h \\ E 



k=l 



+ E J [p(x)- Ph ][V.u(x)- (V. Uh )\ ns ]dx 

k=l fif k 



(64) 
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For the MINI method, its energy error norm which is evaluated on N e trian- 
gles T e G Th is written as 



1/2 



2/i Yl rneas(Q s k )[e(u) - e{u h )] T D ^[e{u) - e{u h )] 



e=l 



T e eT h 



U - Uh\\E 



= < 



> 




(65) 



6 Numerical results 



In this section, we show some numerical results to estimate the efficiency and 
accuracy of the bES-FEM method. For this objective, we use the following bench- 
mark problems to compare between the bES-FEM method and the above methods: 

The first benchmark problem is considered to be the Cook's membrane problem. 
This well-known problem is often used to test for bending behavior and volumetric 
locking 0g, [SH], H3J. Let Q is the convex hull 



The domain f2 is the tapered panel (see Figure 9) whose left boundary is claimed 
, and whose right one is subjected to an in-plane shearing load of IN along the 
y-direction. The two material coefficients are the Young modulus E = 1 and the 
Poisson ratio v = 0.4999999, which follows the nearly incompressible case. For 
this problem, the vertical displacement at the center tip section is found to be 



U = cora;{(0,0),(48,44),(48,60),(0,44)}. 



18.5002 [38J. 



44 



If, 




IN 



Figure 9. Cook's membrane problem on the primal triangulation. 
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In Figures 10, 11 and 12, we represent the comparison of numerical results between 
the bES-FEM and other methods for the Cook's membrane problem. 
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Figure 10. Convergence of displacement tip for Cook's membrane problem (v = 0.4999999), 
where the reference solution is denoted by ref (the solid black line). 

As shown in Figure 10, it is evident that the bES-FEM method is the most accurate 
than the other methods. Moreover, both ES-FEM and Q4 methods are subjected 
to volumetric locking. The NS-FEM method gives an upper bound solution and 
procedures quite well displacement tip without oscillations. Unfortunately, this 
method can no longer be guaranteed the stability of the pressure (see Figure 11). 
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Figure 11. Pressure along the line y (x = 24) for Cook's membrane problem (v = 0.4999999). 

Figure 11 illustrates the pressure distributions whose both the MINI element and 
the bES-FEM method are stable, and whose the NS-FEM method exhibits oscilla- 
tions. However, for instability of the NS-FEM method, it is worth to note that this 
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problem is significantly improved (see Figure 12), when we use the finer triangular 
mesh 16 x 16. Besides, the two figures 11 and 12 indicate that the pressure results 
of both the ES-FEM method is very unstable. 
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Figure 12. Pressure along the line y (x = 24) for Cook's membrane problem (v = 0. 



The next benchmark problem, considered in |10j and [46], is a cylindrical pipe 
subjected to an inner pressure p = 8kN/m 2 , where its internal radius and external 
radius are a = lm b = 2m, respectively (see Figure 13). 





xxn 



X 



Figure 13. A cylindrical pipe subjected to an inner pressure and its quarter model 
with symmetric conditions imposed on the left and bottom edges. 

Due to the axisymmetric characteristic of the problem, we model the upper right 
quadrant of the pipe imposed symmetric conditions on the left, bottom edges, 
and a free traction on its outer boundary. This problem is interested in the nearly 
incompressible case, i.e the Poisson ratio v is close to 0.5, to evaluate frame invari- 
ance and plane bending. Its domain is meshed by 3-node triangular and 4-node 
quadrilateral elements as shown in Figure 14. 
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a) 16x8x2 triangular elements b) 16x8 quadrilateral elements 
Figure 14. a) domain discretization of 3-node triangular b) 4-node quadrilateral elements. 

Additionally, the cylindrical pipe problem, which satisfies plane strain conditions 
and Young's modulus E = 2 1000A; iV/m, has its exact solution for the radial and 
tangential exact displacement [53] 



u r (r) 



(1 + v)a 2 p 
E(b 2 - a 2 ) 



(1 - 2v) + 
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and for the radial and tangential exact stresses 



and u, 







(66) 



a r (r) 



2 

a p 



b 2 



b 2 



2 

a p 



b 2 



1 + 



a 



rip 



(67) 



In two equations (66) and (67), (r, if) are the polar coordinates, and (p is measured 
counter-clockwise from the positive x-axis. 



The MINI, NS-FEM, bES-FEM methods are used to investigate the convergence 
rate of this problem with the Poisson ration v = 0.4999999 whose results illustrate 
in Figure 15, 16 and 17 



-4 r 
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Figure 15. Displacement error norms of the bES-FEM method in comparison with the NS-FEM 
and MINI method for the cylindrical pipe under nearly incompressible condition {y = 0.4999999). 
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Figure 16. Pressure error norms of the bES-FEM method in comparison with the NS-FEM 

and MINI method for the cylindrical pipe under nearly incompressible condition (y = 0.4999999). 
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Figure 17. Energy error norms of the bES-FEM method in comparison with the NS-FEM 

and MINI method for the cylindrical pipe under nearly incompressible condition (y = 0.4999999). 

According to Figures 15 and 16, the two convergence rates in both the displacement 
and the pressure error norms of the bES-FEM method are very high (> 1.95). With 
the MINI and NS-FEM methods, their convergence rates in the displacement error 
norm are close to 2, but their convergence rates in the pressure error norm are not 
as high as that of the bES-FEM method. Moreover, both the displacement and 
the pressure error norms of the bES-FEM method are more accuracy than ones 
of the MINI and NS-FEM methods. Besides, Figure 17 indicates that the error 
in energy norm for the bES-FEM method is smaller than those of the MINI and 
NS-FEM methods. 
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7 Conclusion 



In this work, we have presented an edge-based smoothed element method enriched 
with bubble functions (bES-FEM) for nearly incompressible 2D-3D elasticity. The 
method help soften the bilinear form allowing the weakened weak (W2) procedure 
to get a quality solution. For the bES-FEM method, we show that the uniform 
inf-suf condition satisfies. Besides, the degree of freedom for the pressure can be 
computed by linear combinations depending on the displacement, then the mixed 
displacement-pressure formulation is based only on displacements. 
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